\documentclass[a4paper]{article}
\usepackage{amsmath,amssymb,amsfonts,latexsym,graphicx}
\usepackage[boldfont,slantfont]{xeCJK}
\usepackage{xcolor}
\usepackage{enumerate}
\setCJKmainfont{AR PL UMing CN}


\parindent=0mm
\textheight=20cm \textwidth=14cm
\parskip=0mm
\renewcommand\baselinestretch{1.2}
\oddsidemargin30pt \evensidemargin30pt


\begin{document}

\baselineskip 18pt
\parskip 10pt
\parindent 0em

\def\chntoday{\the\year~年~\the\month~月~\the\day~日}

\renewcommand\figurename{图}\renewcommand\tablename{表}
\renewcommand{\refname}{参考文献}
\renewcommand\abstractname{\bf\large 摘~~要}

\title{计算流体力学课程上机作业(三)}
\author{叶时炜}
\date{\chntoday}
\maketitle

\begin{abstract}
用Harten 的TVD格式计算前台阶问题。
\end{abstract}

\section{问题}

已知二维Euler方程组：
\begin{displaymath}\label{eq:euler}
  \left\{
  \begin{array}{c}
    \left(
    \begin{array}{c}
      \rho\\
      \rho u\\
      \rho v\\
      E\\
    \end{array}
    \right)_t + 
    \left(
    \begin{array}{c}
      \rho u\\
      \rho u^2 +p\\
      \rho uv\\
      u(E+p)\\
    \end{array}
    \right)_x 
    +\left(
    \begin{array}{c}
      \rho v\\
      \rho uv\\
      \rho v^2 +p\\
      v(E+p)\\
    \end{array}
    \right)_y=0,\\
    p=(\gamma-1)(E-\frac{1}{2}\rho u^2-\frac{1}{2}\rho v^2),\gamma=1.4\\
  \end{array}
  \right.
\end{displaymath}
\[
U =  
    \left(
    \begin{array}{c}
      \rho\\
      \rho u\\
      \rho v\\
      E\\
    \end{array}
    \right)
F=    \left(
    \begin{array}{c}
      \rho u\\
      \rho u^2 +p\\
      \rho uv\\
      u(E+p)\\
    \end{array}
    \right),
    G=\left(
    \begin{array}{c}
      \rho V\\
      \rho uv\\
      \rho v^2 +p\\
      v(E+p)\\
    \end{array}
    \right)
\]
(由Mathematica符号计算得出,和课件上面的不太一样。)：
\[A=\frac{\partial F}{\partial U}=
\left(
\begin{array}{cccc}
 0 & 1 & 0 & 0 \\
 \frac{1}{2} \left(u^2 (-3+\gamma )+v^2 (-1+\gamma )\right) & -u (-3+\gamma ) & v-v \gamma  & -1+\gamma  \\
 -u v & v & u & 0 \\
u\left(\frac{1}{2}\text{  }\left(u^2+v^2\right) (-1+\gamma )+H\right) &H-u^2 (-1+\gamma ) & -u v (-1+\gamma ) & u \gamma 
\end{array}
\right)
H = \frac{E+p}{\rho}
\]
\[B=\frac{\partial G}{\partial U}=
\left(
\begin{array}{cccc}
 0 & 0 & 1 & 0 \\
 -u v & v & u & 0 \\
 \frac{1}{2} \left(v^2 (-3+\gamma )+u^2 (-1+\gamma )\right) & u-u \gamma  & -v (-3+\gamma ) & -1+\gamma  \\
 v\left(\frac{1}{2}\text{  }\left(u^2+v^2\right) (-1+\gamma )+H\right) & -u v (-1+\gamma ) & -v^2 (-1+\gamma )+H & v \gamma 
\end{array}
\right)
H = \frac{E+p}{\rho}
\]
A的特征值与特征向量(由Mathematica符号计算得出，并通过适当变换使分母不出现$u,v$)：
\[\lambda_1=u-a,\lambda_2=\lambda_3=u,\lambda_4=u+a,a=\sqrt{\frac{\gamma p}{\rho}}\]
\[ R^{(1)}=\left( \begin{array}{c}
1\\
u-a\\
v\\
H-au\\
\end{array}\right)，
 R^{(2)}=\left( 
\begin{array}{c}
 0 \\
 0 \\
 1 \\
 v
\end{array}
\right)，
 R^{(3)}=\left(
\begin{array}{c}
 1 \\
 u \\
 v\\
 \frac{u^2+v^2}{2} 
\end{array}
\right)，
 R^{(4)}=\left( \begin{array}{c}
1\\
u+a\\
v\\
H+au\\
\end{array}\right)。
\]
\[L=\left(
\begin{array}{cccc}
 \frac{2 a u+\left(u^2+v^2\right) (-1+\gamma )}{4 a^2} & -\frac{a+u (-1+\gamma )}{2 a^2} & -\frac{v (-1+\gamma )}{2 a^2} & \frac{-1+\gamma }{2 a^2} \\
 -v & 0 & 1 & 0 \\
 \frac{2 a^2-\left(u^2+v^2\right) (-1+\gamma )}{2 a^2} & \frac{u (-1+\gamma )}{a^2} & \frac{v (-1+\gamma )}{a^2} & \frac{1-\gamma }{a^2} \\
 \frac{-2 a u+\left(u^2+v^2\right) (-1+\gamma )}{4 a^2} & \frac{a+u-u \gamma }{2 a^2} & -\frac{v (-1+\gamma )}{2 a^2} & \frac{-1+\gamma }{2 a^2}
\end{array}
\right)\]
B的特征值与特征向量(由Mathematica符号计算得出，并通过适当变换使分母不出现$u,v$)：
\[\lambda_1=v-a,\lambda_2=\lambda_3=v,\lambda_4=v+a,a=\sqrt{\frac{\gamma p}{\rho}}\]
\[ R^{(1)}=\left( \begin{array}{c}
1\\
u\\
v-a\\
H-av\\
\end{array}\right)，
 R^{(2)}=\left( 
\begin{array}{c}
 0 \\
 1 \\
 0 \\
 u
\end{array}
\right)，
 R^{(3)}=\left(
\begin{array}{c}
 1 \\
 u\\
 v \\
 \frac{u^2+v^2}{2 } 
\end{array}
\right)，
 R^{(4)}=\left( \begin{array}{c}
1\\
u\\
v+a\\
H+av\\
\end{array}\right)。
\]
\[
L=\left(
\begin{array}{cccc}
 \frac{v (2 a+v (-1+\gamma ))+u^2 (-1+\gamma )}{4 a^2} & -\frac{u (-1+\gamma )}{2 a^2} & -\frac{a+v (-1+\gamma )}{2 a^2} & \frac{-1+\gamma }{2 a^2} \\
 -u & 1 & 0 & 0 \\
 \frac{2 a^2-\left(u^2+v^2\right) (-1+\gamma )}{2 a^2} & \frac{u (-1+\gamma )}{a^2} & \frac{v (-1+\gamma )}{a^2} & \frac{1-\gamma }{a^2} \\
 \frac{v (-2 a+v (-1+\gamma ))+u^2 (-1+\gamma )}{4 a^2} & -\frac{u (-1+\gamma )}{2 a^2} & \frac{a+v-v \gamma }{2 a^2} & \frac{-1+\gamma }{2 a^2}
\end{array}
\right)
\]
前台阶问题，几何形状如图\ref{pic:step}。初始时刻，区域内充满均匀流$(\rho,u,v,p)=(1.4,3,0,1)$,
输出时刻为$t=4$.上下边界为反射边界，左边界为入流边界；右边界为出流。
\begin{figure}\centering
  \includegraphics[width=14cm]{./step.png}
    \caption{前台阶问题的几何形状}\label{pic:step}
\end{figure}
\section{数值方法}
\subsection{Spang-type dimensional splitting}
\begin{equation}
  v^{n+2}=Lv^n
\end{equation}
\begin{equation}
  L = L_xL_yL_yL_x
\end{equation}
其中$L_x$, $L_y$为下列一维方程的格式：
\begin{equation}
  L_x:w_t+f(w)_x=0,  L_y:w_t+g(w)_y=0,
\end{equation}
如果$L_x$和$L_y$是上面方程的二阶近似，那么$L$就是下面方程的二阶近似：
\begin{equation}
  w_t+f(w)_x+g(w)_y=0,
\end{equation}


\subsection{Harten 的TVD格式}
ULT1C格式\cite{thz03}。
\begin{equation}
v^{n+1}=v^{n}-\lambda(\bar{f}_{j+1/2}-\bar{f}_{j-1/2}),
\end{equation}
其中：
\begin{align}
  \bar{f}_{j+1/2}=&\frac{1}{2}[f(v_j)+f(v_{j+1})-\frac{1}{\lambda}\sum_{k=1}^{m}\beta_{j+1/2}^{k}R_{j+1/2}^{k}]\\
  \beta_{j+1/2}^k=&Q(\nu_{j+1/2}^k+(1+2\hat{\theta}_{j+1/2})\gamma_{j+1/2}^k)\alpha_{j+1/2}^k
  -(1+2\hat{\theta}_{j+1/2})(g_j^k+g_{j+1}^k)\\
  \nu_{j+1/2}=&\lambda a^k{v_{j+1/2}},a\,is\,eigen\,value\\
  g_{j}^k=&s_{j+1/2}^kmax(0,min(|\tilde{g}_{j+1/2}^k|,\tilde{g}_{j-1/2}^ks_{j+1/2}^k))\\
  s_{j+1/2}^k=&sign(\tilde{g}_{j+1/2}^k)\\
  \tilde{g}_{j+1/2}^k=&\frac{1}{2}(Q(\nu_{j+1/2}^k)-(\nu_{j+1/2}^k)^2)\alpha_{j+1/2}^k\\
  \gamma_{j+1/2}^k=&\left\{
  \begin{array}{l l}
    (g_{j+1}^k-g_{j}^k)/\alpha_{j+1/2}^k, & when\, \alpha_{j+1/2}^k \neq 0\\
    0,&when\, \alpha_{j+1/2}^k = 0\\
    \end{array}
  \right.\\
  Q(x)=&x^2+\frac{1}{4}, for\,all\,k\\
  \hat{\theta}_{j+1/2}^k=&max(\theta_{j}^k,\theta_{j+1}^k), for\,all\,k\\
  \theta_{j}^k=&|\alpha_{j+1/2}^k-\alpha_{i-1/2}^k|/(|\alpha_{j+1/2}^k|+|\alpha_{i-1/2}^k|)
\end{align}
\begin{align}
  v_{j+1/2}=&\frac{v_{j}+v_{j+1}}{2}\\
  \alpha_{j+1/2}=&L_{j+1/2}(v_{j+1}-v_j)\\
\end{align}
$L_{j+1/2},R_{j+1/2}$是系数矩阵的左右特征向量矩阵。方程中若直接取$\hat{\theta}=0$就是ULT\cite{thz03}格式。

\subsection{边界条件}
建立矩形网格使边界正好在两个网格节点中间，迭代需要额外两个影子节点。


边界条件\cite{thz04}：
\begin{description}
  \item[入流边界]将两个影子节点设为入流初始值。
  \item[反射边界]将影子节点的与反射面法向的速度取作镜像位置的点的相反值。
  \item[出流边界]直接将影子节点的值取作边界节点的值。
\end{description}
台阶角特殊处理\cite{thz04}。

\section{数值实验}
我用matlab实现了上述算法，具体版本为2010B，linux。
\subsection{程序简述}
\begin{enumerate}
  \item LTVD.m: LTVD(Uin,boundary,lambda,method)一维TVD格式的迭代步，根据参数输入，确定。
    \begin{description}
      \item[Uin]前一步数据。
      \item[boundary] 边界。
      \item[lambda] $\lambda = \frac{\tau}{dx}$
      \item[method] 方法ULT1 or ULT1C，ULT1C没有调好。
    \end{description}
  \item solvestep.m: solvestep(CFL,n,opt,method)
    \begin{description}
      \item[CFL]CFL系数。
      \item[n] 纵向网格数，必须是5的倍数。
      \item[opt] opt = 1时，动态画出计算过程中的$\rho$，会浪费很多计算时间。
      \item[method] 方法ULT1 or ULT1C
    \end{description}


    sovlestep调用LTVD.m计算台阶问题,计算过程中会输出计算进度，按$\frac{t}{T}$计算,$t$为当前计算时刻，$T$为最后显示时刻。
  \item main.m: 此程序调用solvestep.m计算不同网格大小的结果。
\end{enumerate}

\subsection{实验结果}
网格数取不同值得到以下结果。
\begin{enumerate}
  \item n=10 ，如图\ref{pic10}
  \item n=20 ，如图\ref{pic20}
  \item n=40 ，如图\ref{pic40}
\end{enumerate}


计算时间(以下为main函数的输出结果)：
\begin{verbatim}
n=10, starting...
...................................................100.0%

 n=10, Done. Spends 45.1678 seconds
n=20, starting...
...................................................100.0%

 n=20, Done. Spends 290.5631 seconds
n=40, starting...
...................................................100.0%

 n=40, Done. Spends 2144.2951 seconds
\end{verbatim}
\begin{figure}[htbp]
  \centering
\includegraphics[width=14cm]{../n=10-a.eps}
\caption{n=10时的计算结果}
\label{pic10}
\end{figure}
\begin{figure}[htbp]
  \centering
\includegraphics[width=14cm]{../n=20.eps}
\caption{n=20时的计算结果}
\label{pic20}
\end{figure}
\begin{figure}[htbp]
  \centering
\includegraphics[width=14cm]{../n=40.eps}
\caption{n=40时的计算结果}
\label{pic40}
\end{figure}


\begin{thebibliography}{99}
% adaptive
\bibitem{thz01}
李治平，《偏微分方程数值解讲义》，2009
\bibitem{thz02}
汤华中，《计算流体力学讲义》，2011
\bibitem{thz03}
Ami Harten, High Resolution Schemes for Hyperbolic Conservation Laws, JOURNAL OF COMPUTATIONAL PHYSICS 135, 260–278 (1997)
\bibitem{thz04}
Richard Sanders, High Resolution Staggered Mesh Approach for Nonlinear
Hyperbolic Systems of Conservation Laws, JOURNAL OF COMPUTATIONAL PHYSICS 101, 314329 (1992)
\end{thebibliography}
\end{document}
